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ABSTRACT 

Context. The solution of the nonlocal thermodynamical equilibrium (non-LTE) radiative transfer equation usually relies on stationary 
iterative methods, which may falsely converge in some cases. Furthermore, these methods are often unable to handle large-scale 
systems, such as molecular spectra emerging from, for example, cool stellar atmospheres. 

Aims. Our objective is to develop a new method, which aims to circumvent these problems, using nonstationary numerical techniques 
and taking advantage of parallel computers. 

Methods. The technique we develop may be seen as a generalization of the coupled escape probability method. It solves the statistical 
equilibrium equations in all layers of a discretized model simultaneously. The numerical scheme adopted is based on the generalized 
minimum residual method. 

Results. The code has already been applied to the special case of the water spectrum in a red supergiant stellar atmosphere. This 
demonstrates the fast convergence of this method, and opens the way to a wide variety of astrophysical problems. 
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1. Introduction 

Radiative transfer plays a central role in astrophysics. The solu¬ 
tion of the radiative transfer equation (hereafter RTE) is required 
to properly interpret any spectrum observed from an astrophys¬ 
ical object. Eurthermore the radiation field often represents an 
important ingredient in hydrodynamics, for instance, through the 
energy budget (radiative heating and/or cooling) and the radia¬ 
tive pressure. However, as we discuss below, while the standard 
expression of the RTE is an ordinary differential equation of the 
1st order, it may be highly nonlinear and nonlocal, especially 
because of scattering processes. The solution of the RTE thus re¬ 
quires specific procedures. Eurthermore, with the development 
of infrared (IR) astronomy (e.g., from CRIRES@VLT in the 
near-IR to ALMA through HIEI/Herschel), one is now often con¬ 
fronted with very large-scale problems resulting from molecular 
spectra, which most techniques developed so far cannot handle. 

The techniques for solving the RTE may be divided into a 
few strategies. Among approximate methods, one that is widely 


used, relies on the escape probability technique (e.g.. Castor 
|1970^ , which does not take the nonlocality of the RTE into ac¬ 
count. Secondly, Monte Carlo methods can be used because of 
the probabilistic nature of scattering processes. These methods 
can successfully be adapted to any geometry, but may be un¬ 
adapted when large optical depths are met (see, e.g., lJuvela 


|2005| for a recent review). We do not discuss these two classes 
of methods hereafter. 


To properly treat the nonlocality and the nonlinearity of the 
RTE, it is usually rewritten in the form of a diffusion integral 
equation, which is then solved iterativeljj^ This integral form is 
usually referred to as the A operator. The form involves a ker¬ 
nel, which determines the mathematical nature of the integral 
equation. Eor example, in the case of an assumed plane parallel 
geometry, the kernel involves the exponential integral function, 
and the RTE then becomes a weakly singular Eredholm equation 
of the second type (this singularity is an important aspect, which 
will be treated below). It can then be formally solved thanks to 
appropriate algorithms, such as the Atkinson algorithm ( |Ahues 
|et al.|20()2) l. Since the inversion of the A operator is usually nu¬ 
merically prohibitive, one needs an iteration scheme. It is well 
known that the iteration with the complete A operator at best 
converges very slowly or even stabilizes far from the real solu¬ 
tion. A major step in solving the RTE has been the introduction 
of the so-called approximate lambda iteration (ALI), initially by 
|Cannon| ( |1973| l, then called the operator splitting method. The 
ALI is in fact an application of the Jacobi method to the RTE, as 
shown by Olson et al. ( |1986| l, who greatly improved the math¬ 
ematical understanding of the problem. In these methods, the 
choice of the approximate operator (denoted A*) is crucial. The 


' We already emphasize here that solving either the RTE explicitly or 
the statistical equilibrium equation is strictly identical, and both equa¬ 
tions are of similar mathematical nature, as shown, for instance, in the 
Appendix of Hauschildt et al. jl995| 
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most common choice for A* is the diagonal part of A, which 
ensures a rapid convergence if A is diagonally dominant matrix 
(Gershgorin’s Theorem). ALI has been applied to a large vari- 


ety of problems, including 3D polarized transfer ( Stepan & Tru 
jillo Bueno|2013 1. This latter work also emphasizes the potential 


progress that can be made by taking advantage of massively par¬ 
allel computing. 

Meanwhile, new mathematical techniques have been applied 
to the RTE problem, especially the Gauss-Seidel method and the 
successive over relaxation (SOR) method, which significantly 
improve the convergence properties ( [Trujillo Bueno & Fabiani| 
|Bendicho |1995| l that are determined by the spectral radius of 
their amplification matrix. Besides these, complete linearization 
methods have also been developed, the well-known case being 
MULTI ( [Scharmer & Carlsson||1985] l. Despite apparent differ¬ 
ent formulations, the two approaches are essentially the same, as 
shown by [Socas-Navarro & Trujillo Bueno| ( |1997| l. 

Indeed, all these methods are based on linearization tech¬ 
niques. In particular, the dependence of the radiation field on the 
populations of the energy levels is approximated. Furthermore, 
all these methods are stationary (in the case of SOR, this is true 
if the relaxation parameter is constant over iterations), meaning 
that the spectral radius is constant. As the error vector e™ at the 
m''^ iteration is proportional to the power m of the spectral radius 
of the amplification matrix Idmaxl'", the convergence may thus 
be very slow in some cases. For the Lambda iteration, the spec¬ 
tral radius is close to 1 when the photon destruction probability 
e 1 and/or the total optical depth is » 1, explaining why it 
falsely converges. This is illustrated in Figs.[2and[^ where the 
classical Eddington problem (two-level atom with constant e) 
is solved with A iteration (LI) and ALI, respectively. The false 
convergence of LI and the very low convergence of ALI are illus¬ 
trated in cases B. Indeed, their rate of convergence is very slow 
for strongly nonlocal thermodynamical equilibrium (non-LTE) 
problems (non-LTE parameter e « 1, e being the fraction of 
the source function given by the Planck function, or the photon 
destruction probability) and to the number of points per decade 
of optical depth. In the latter case (i.e., a fine sampling of the 
stellar atmosphere), this problem of slow convergence is well 
explained by |Fabiani Bendicho et aL| ( |1997| l, who also showed 
that the multigrid method can nicely circumvent this problem. 
This fine sampling is however rarely required, especially as clas¬ 
sical model atmospheres are static. An excellent discussion of 
the general problem of false convergence is given in Chapter 13 
Hubeny & Mihalas ( 2014| ). 


in 


The iterative schemes are usually combined with a technique 
of acceleration of convergence ( Auer|199T] l. These fall into two 
main categories. The first is based on the minimization of the 
residual, which is the case of the Ng acceleration ( |Ng 1974| l. 
The second class relies on minimization with respect to a set of 
conjugate vectors. This was first applied to the RTE by ( |Klein| 
|et al.|198^ ; see also |Dickel & Auer|1994[ ). 

The second class of techniques (the minimization with re¬ 
spect to a set of conjugate vectors) leads naturally to more gen¬ 
eral, nonlinear, nonstationary methods, such as the conjugate 
gradient and the generalized minimum residual method (GM- 
RES). To our knowledge, only the conjugate gradient has been 
applied to the RTE, by Paletou & Anterrieu ( |2009 1 . 

The aim of our work is to develop a new nonstationary 
method for solving the RTE, which fully takes its nonlinearity 
into account. This new method thus opens the possibility to ex¬ 
plore the mathematical consequences of this intrinsic property 
of the RTE. Furthermore, circumventing the possible slow or 
false convergence of stationary methods should allow an appro- 


Article number, page 2 of 13 


Spectral Radii (P ]forA —0) Two levels benchmark 

- 0.00 ^ 

-0.05 S 
of 

—0.10 be 

_o 

-0.15 


- 0.00 
-0.50 g 
- 1.00 5 
-1.50 S 
- 2.00 

0 0.2 0.4 0.6 0.8 1 10-‘ 10" 10‘ 10^ 10" 

Destruction probability s Optical Depth fv 



Fig. 1. Spectral radius of amplification matrix for the A iteration 
method. Convergence is shown for two spectral radii, marked with let¬ 
ters A (upper panel) and B (lower panel) in the right-hand panels. The 
thin black lines correspond to successive iterations; the true solution is 
given by the orange line. 


Spectral Radii (P | for A — A(t)I) Two levels benchmark 


Gh 

c 

'o 

cu 


512 

128 

64 

32 

16 

8 

4 

2 



1 



- 0.00 

-0.05 5a 

- 0.10 ^ 
_o 

-0.15 ^ 

- 0.00 
-0.50 ^ 
- 1.00 5 
-1.50 .o' 
- 2.00 


Destruction probability £ 


Optical Depth Ty 


Fig. 2. Same as Figure[^ except for ALI with A* = diagonal of A. 


priate treatment of some specific cases, such as, e.g., shocked 
regions, which require a fine sampling and present a large dy¬ 
namics of optical depth values. Finally, as the coding strategy of 
this method has been conceived as a parallel code from the be¬ 
ginning, it is applicable to large-scale systems such as molecular 
spectra, and this was indeed our initial motivation (a posteriori 
parallelization is often less efficient, or even impossible). 

The derivation of the equations, including the Jacobian ma¬ 
trix coefficients, used in the nonlinear method is presented in 
section 2. The implementation of this method and its paralleliza¬ 
tion are described in section 3. The application of the code to a 
specific problem facilitates the determination of its convergence 
properties, as shown in section 4. 


2. Global strategy 

Contrary to the two-level atom problem with a constant destruc¬ 
tion probability, which can be directly solved by direct inversion 
(or through the analytic Eddington solution). Multi-level atom 
or molecule problems require an iterative determination for two 
reasons. Firstly, the multi-level problem is nonlinear. Secondly, 
the radiative transfer equation is nonlocal, and one has to prop¬ 
agate the modification of the field during iterations. As we dis¬ 
cussed in the previous section, however, this stationary iteration 
may be problematic for large-scale systems. The basic concept 
of the strategy we adopt is to suppress the iterations due to the 
nonlocality of the RTE by considering the full system (spatially) 
and computing the radiation field exactly (assuming given popu¬ 
lations) with an analytic formulation of the mean intensity field, 
which explicitly depends on populations. 
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2.1. Formulation of the problem 


For clarity, we hereafter consider the special case of plane paral 
lei geometry. In this context the ID RTF is 

dly 




ds 
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( 1 ) 


where ly is the specific intensity, rj andx are the emission and ex¬ 
tinction coefficients, and s is the geometrical path. We explicitly 
separate the continuum and the line processes (all demonstra¬ 
tions and symbols can be found in Appendix and respec¬ 
tively). Assuming complete redistribution, we introduce classi¬ 
cal expressions for line processes. 
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where u and I refer to the upper and lower levels, respectively. 
Hereafter, if the distinction between upper and lower levels is 
not required, the indices i and j are used. The normalized line 
profile is and A and B are the Einstein coefficients. 

Following the strategy of Gonzalez Garcia et al. ( 2008| l, we 
then introduce a unique optical depth scale for all frequencies, 
e.g,. the commonly used ry scale (optical depth at 500 nm). Eq. 
0 then becomes 

dL 
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The notation T means that the quantity is expressed per hydro¬ 
gen atom. Details can be found in Appendix [A| As implicitly as¬ 
sumed in Eq.|3 we do not consider line overlap. This assumption 
is justified in the case of infrared molecular lines. Otherwise, a 
summation over transitions (at least adjacent ones) would be re¬ 
quired. 

This formulation depends explicitly on the populations of 
each energy level of the species considered. These populations 
are solutions of the statistical equilibrium equation system 
which is for each level n,, the detailed balance of incoming anc 
outgoing (de-)excitation processes, defined by radiative (Ein¬ 
stein coefficients), and collisional rates (ticoi being the local den¬ 
sity of the colliding partner). Thus the statistical equilibrium 
equations take the form 


drii 

dt 




j<‘ 




^ ^ rij (BjiSf ji + ricoiCj^ 


j>‘ 


j*i 


= 0, 


(5) 


^ At steady-state. Chemical formation and destruction terms are not 
included. 


where the mean radiation field J/j/Ty) = Sf jAtv) depends on the 
solution of the RTE, 

2 ^OO 

ipy(Tv)j Iy(Ty,p)dpdv. (6) 

Because this system of equations is degenerate, one of them is 
replaced by the normalization equation where is 

the total density of the species k. 

To express each equation of the statistical equilibrium sys¬ 
tem with the same unknowns as in Eq. we rewrite the statisti¬ 
cal equilibrium with the variables v, defined in Eq. Moreover, 
for numerical reasons, we make a variable change for the mean 
intensity. 


2hv^. 

JijiTv) = -^ZijiTv). (7) 

This leads to a matrix formulation of the statistical equilib¬ 
rium, with an upper triangular part corresponding to incoming 
processes, and a lower triangular part corresponding to outgo¬ 
ing processes. Finally, thanks to the Einstein relations {B^i = 
(c^l(2hvl,))Aui and B/„ = (gulgi)Bui), the system of equations 
0 becomes 


^ gi^'ii} 2 ^ 0 ') ^ gj^i'^ij ^ ricoigiC, 

j<i j>i i*i 

^ ^jgj-^ij^ij ^ ^jgj^ji ^Ji) ^ ricolgiC, 
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= 0 . 

( 8 ) 


Because radiative transfer is a nonlocal problem, an integral for¬ 
mulation of the mean radiation field with an explicit dependence 
on populations is required. 


2.2. The mean radiation field 

To derive the mean radiation field, the first step consists of deriv¬ 
ing the formal solution of the RTE. A classical way to obtain an 
integral form of the formal solution of a differential equation is to 
use the Green’s functions This method consists of replacing 
the source term by a Dirac function. Taking the boundary condi¬ 
tions into account, the equation is then solved using its Laplace 
transform. The complete solution is then obtained after integra¬ 
tion of this Green’s function, that is, performing a summation 
over each point source, according to the superposition principle 
(which is valid if the differential operator is linear). Thus, the 
integral form of /y(Tv,ju) can be obtained by the determination 
of the Green’s function of Eq. 0. This Green’s function is eas¬ 
ily computable for spherical and other geometries using some 
symbolic tools such as Mathematica, making it possible to gen¬ 
eralize the method for 2D or 3D problems. Knowing this Green’s 
function, we can directly obtain the formal integral form of the 
specific intensity, which is 


/y(TV,jU) = f \DijXi{t)^{t)(l).^{t) + ^y{t)Sy{t)\%y{TY\t)dt. 

I -Min -* 

(9) 

We assume the line profile does not depend on p. This expres¬ 
sion is thus not applicable in the case of anisotropic scattering 
or in the presence of radiation fields, but is valid in the case of 
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classical static model atmospheres. The formal expression of the 
mean intensity, requires an integration over solid angles 

and frequencies according to Eq We assume that DjjXi^ and 
^ySv are angular independent (physically this corresponds to the 
assumption of isotropic scattering) and frequency independent 
over the line width {^ySy This is a very weak as¬ 

sumption because the continuum is nearly constant on the scale 
of a line width. Then, we can extract the source terms from the 
integral over solid angles and frequencies. Moreover, we sepa¬ 
rate line and continuum contributions by splitting the integral, 
which leads to 
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( 10 ) 


Then, after determination of the Green’s function in a given ge¬ 
ometry, one obtains an integral form for the different contributors 
to the radiation field. We assume that the continuum contribution 
is in LTE (Sy.j = By..(T)). Moreover, we consider, as boundary 
conditions, an incoming field 7®’“ from the inner atmosphere. 
Then, 


jijiry) = J5’“(tv) + + J/^™(tv) 

where. 
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^ij{t)By..(T (t))Li [Xi, Xj\(t, TY)dt (12b) 


Xi{t)((f)K\ [Xi, Xj\{t, Ty)dt, (12c) 


where /3_[x,, Xj](Tv) is the escape probability from position Ty 
in the atmosphere, 

/ 3 -[Xi,Xj](Ty) = - j (p'JiTy) XE2(y°'^[Xi,Xj](Ty))dv. (13) 

The kernels Li [xj, Xj]{t, Ty) and [x,-, Xj]{t, Ty) are, respec¬ 
tively: 

1 i ■ 

Ll[X/,Xj](f,Tv) = - I <^v(tv) 

xEi (-s(f, Tv)(T'^°‘[xi, Xj](t) - x^](tv))) dv (14a) 


Ki [xi,Xi](t, 
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where the total optical depth is given by 


(14b) 
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) dt 
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(^y(t) + Eu^(t)(p‘y^(t)[xj(t) - Xi(0])dt. 
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We do not include an external radiation field. In the plane par¬ 
allel geometry considered here, it would require us to take ad¬ 
ditional effects such as limb darkening in realistic cases into ac¬ 
count (e.g., illumination by a planet or a companion). We have 
introduced the function e(t, Ty) = 1 if f < Ty and - 1 if f > Ty 
to split the integration over p according to its sign. The physical 
meaning of the kernels L\{t, Ty) and Ki{t, Ty) can be understood 
as the probability density for photons, emitted in a layer f, to 
contribute to the line excitation in the layer Ty. 

We now have an analytic expression of the mean intensity with 
an explicit dependence on populations. We now have to include 
this expression in the equations of statistical equilibrium to close 
the system and form a system of coupled integral equations (a 
Eredholm-like equation system). 

2.3. The multi-zone statistical equilibrium equation system 

The inversion of the system of statistical equilibrium equations 
allows us to obtain the populations also including the formal so¬ 
lution of RTE. But because it is nonlocal, the knowledge of S 
requires the knowledge of the populations everywhere in the at¬ 
mosphere because of the integral over the optical depth t in Eq. 

GH). 

Eor this reason, one needs to solve the statistical equilibrium 
everywhere in the atmosphere simultaneously. The first step con¬ 
sists of rewriting Eq. 0 in a matrix form. A vector x(Ty) con¬ 
taining all the populations of the considered species in a given 
layer, normalized to the density in that layer, is built. The size of 
this vector is thus equal to N, the number of energy levels of the 
considered species. We introduce the rate matrix 

R[x](Ty)=A^o(Ll+Z[x](Ty)f 

-I- diag(g) ■ [a o Z[x](Ty)] ■ diag(g)^' -H C^(Ty). (16) 

Where Li is a lower triangular matrix of ones, A is a vector 
containing the spontaneous de-excitation coefficients (A„/), 
Z is a matrix containing the radiation field (see Eq. |^, C a 
matrix containing the collisional rates (ticoiCij), and g is a vector 
containing the statistical weights (g,). 

Now the statistical equilibrium system can be written as, 

M^[x](Ty)x(Ty) = 0, (17) 

with 

M[x](Ty) = [diag(u ■ R[x](Ty)) - R[x](Ty)] ■ diag(g), (18) 

where u is a ones vector. After discretization of the medium into 
P layers indexed by /, the vectorial function x(tv) becomes x’ 
and equation ( [T6l l becomes 

R'(x',x2,...,xO = A^o(Li+ZV,xX...,xOf 

-H diag(g) ■ [a o Z'(xXxX ...,x^)] ■ diag(g)^' H- [c']^ (19) 

as the radiation field, and thus the matrix Z depends implicitly 
on the populations in all layers , and the statistical equilibrium 
system can be rewritten, 

[m'(x',xX...,x^)]^x'= 0. (20) 

To take care of the coupling between the layers and close the 
system to solve Eq. (|20ll for all I in a unique step, we build the 
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multi-zone statistical equilibrium system: 
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( 21 ) 


We introduce the vector X = (x’,x^, and the block- 

diagonal matrix T to rewrite Eq. as a multi-D nonlinear 
function, 


[r(X)]^x = F(X) = 0. 


( 22 ) 


With an initial Xq within the convergence radius, we can solve 
this nonlinear equation with a classical scheme such as the New¬ 
ton method. An important question, however, is to determine 
whether the above system has a solution and whether the so¬ 
lution we expect to get is physical (unique and positive). 


2.4. Existence, uniqueness and positivity of the soiution 

From a physical point of view, one naturally expects a unique 
and positive solution of the RTE. Reflecting on this problem is 
in fact a way to verify whether it is mathematically well formu¬ 
lated. Indeed, as stated above, the coupling between the popula¬ 
tions and the radiation held is nonlinear. Without linearization, 
the solution can be regarded as the equilibrium state of Markov 
processes. Furthermore, because of radiative selection rules and 
collisional propensity rules, many factors are very small, and 
may numerically tend toward zero. The problem is thus poten¬ 
tially ill-conditioned. 

The formulation adopted here is a multi-zone statistical equi¬ 
librium. The problem of uniqueness and positivity of the solu¬ 
tion of the mono-zone statistical equilibrium has been addressed 
by [Damgaard et al. ( |1992| l and |Rybicki| \\991) . As this sys¬ 
tem of equations is degenerate, one equation is usually replaced 
by a normalization condition, namely the conservation equa¬ 
tion, to overcome the singularity of the associated matrix. Then, 
|Damgaard et al.| \\992) show that the solution is positive if all 
the rates are strictly non-zero. Yet certain transition probabili¬ 
ties vanish, physically or numerically. As long as all the levels 
are connected, just one nonzero cofactor guarantees the regular¬ 
ity of the matrix. More generally, Rybicki ( \991\ demonstrated 
that even for the system without a normalization condition, pos¬ 
itive states can link (and be linked) only to positive states. Then, 
if the statistical equilibrium equations can be transformed into 
a set of uncoupled irreducible subproblems with a strictly posi¬ 
tive normalization condition, the solution is positive and unique. 
In other words, the uniqueness and positivity properties of the 
general linear statistical equilibrium equations relate to the un¬ 
derlying connectivity properties of the states. 

In our multi-zone formulation, the coupling between statis¬ 
tical equilibria in every atmospheric layer is strongly nonlinear. 
That is, the system cannot be split into uncoupled sets of equa¬ 
tions, i.e., a set of uncoupled irreducible subproblems. This may 
impact the number of required normalization conditions, while 
only one per atmospheric layer (the conservation equation) is 
available. The uniqueness of the solution may then be demon¬ 
strated through the search for a nonlinear fixed point as the so¬ 
lution of the equilibrium of a Markov chain (R Azerad, priv. 
comm.). A forthcoming paper will be devoted to this. At least, 
as the properties of connectivity still apply, the positivity of the 
solution is guaranteed. 


3. Implementation 

3.1. Noniinear soiving 

In the method developed here, the computation of the statis¬ 
tical equilibrium is based on its exact formulation and explic¬ 
itly includes radiative transfer effects through the mean radiation 
held S in equation ( |2^ . Finding the root of this function natu¬ 
rally leads to the exact solution to the statistical equilibrium, and 
avoids any false convergence that may be obtained with station¬ 
ary iterative methods. 

The computation of the norm of a function, 

F : ^ (23) 

N and P being the number of considered energy levels and the 
number of layers in the model atmosphere, respectively, is equiv¬ 
alent to an optimization problem, that is, the search for the min¬ 
imum of the L 2 norm of the function /, 

min /(X) = iF(X)F^(X) = ^||F(X)||i^. (24) 

The classical way to minimize this function is to move 
along a vector dX in a descending direction, using the Newton’s 
method. Indeed, for bX = X^. - X^_i = ■ F, the gradient 

along the vector bX is 

V/(X,)^bX = [F(X,)^J(X,)] [-J(X,)-'F(X,)] 

= -F(X,)^F(X,) = -||F(X)||i^ < 0. (25) 

Newton’s method is very efficient in approaching the solu¬ 
tion as it has a quadratic convergence. For each successive New¬ 
ton step, we have to solve the linear problem JbXj = -F, where 
k is the index of the Newton iteration. However, with classical 
algebric methods, this operation leads to an O(N^) algorithm, 
which is prohibitive in terms of computational time for large- 
scale systems. A robust and efficient way to accelerate this oper¬ 
ation is to approximate the solution using a faster algorithm. We 
adopt GMRES, which is a subclass of Krylov subspace methods. 

A Krylov subspace method begins with an initial guess bX®. 
At the n* sub-iteration, dX^ is determined through a correction 
in the n* Krylov subspace. 


'K'' = span(R, JR,..., J"''R); R = -F - JbX”. (26) 

The main advantage of this method is that its implementation 
requires only some vector matrix (or vector transposed matrix) 
products. It leads to a complexity in 0{n^), which may be re¬ 
duced to 0{n) if the matrix is sparse. Moreover, the scheme needs 
0{kn) operations, but the algorithm can be restarted to minimize 
this contribution. 

For each Newton step, we thus have to solve the linear 
problem JdX = -F, which will be rewritten in this section as 
Ax = b for simplicity. Stationary iterative methods (Jacobi, 
Gauss-Seidel, SOR) approximate the unknown vector x with a 
scheme x® = + c, where neither B nor c depend upon the 

iteration count k. We use the GMRES method, which is part of a 
class of (nonstationary) methods named Krylov subspace meth¬ 
ods. The GMRES algorithm tries to evaluate the guess x„ e “X), 
(the n-th order Krylov subspace) by solving an optimization 
problem (minimizing the residue ||Ax„ -b||L,) like a least-square 
problem. The n-th Krylov subspace is spanned over the n basis 
vectors, 

'Kn = span(b,Ab,A2b...,A''^'b). (27) 
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Then, the approximate solution is x„ = K„c, where the Krylov 
matrix is K„ = [b,Ab,A2b...,A''-*b] and c is an appropriate 
vector such that 

min ||Ax„ - b||L, = min ||AK„c - b ||/.2 ■ (28) 

This vector is easily computed as it just requires a Matrix-vector 
product. This is particularly well suited for sparse matrices be¬ 
cause this operation evolves roughly as the number of nonzero 
elements. 

A straightforward method to find the least-square solution of 
this problem would be to compute the QR factorization of the 
matrix AK„. Indeed, the normal equation A^Ax - A^b Rx - 
Q^b can be solved by back substitution. However, this is both 
unstable and too expensive (the QR factorization requires 2mn^ 
flops, d - Q^b, 2mn flops, and Rx = d, rP' flops, where m is 
the iterate and n is the order of the Krylov subspace). Moreover 
this natural Krylov basis is ill-conditioned and is thus not a good 
choice. Indeed, the basis vectors have to be linearly independent, 
and the successive multiplications by A lead to vectors, which 
point (fastly) to the dominant eigenvector of A and are thus nu¬ 
merically collinear. For this reason, we build another space , Q„ , 
where the basis vectors are orthogonal (this orthogonalization of 
Krylov bases is named the Arnoldi method). During the Arnoldi 
step to build the orthogonal basis, the matrix H of the orthog¬ 
onalization elements is built. This matrix has a special shape 
(Hessenberg shape), which is close to triangular (one subdiag¬ 
onal more). Our minimization problem can be rewritten as 

min ||Ax„ - bllz., = min ||AQ„y - b||z,j = min ||Q„+iH„+iy - bllz.^, 

(29) 


where y designs the new unknown (Q„y = x„). 

Furthermore, we have to combine this method with a global 
convergence strategy. Indeed, as one gets close to the solution, 
the choice of the descending direction may lead to spurious ef¬ 
fects, especially if dX becomes very large. Then, to ensure the 
global convergence of the scheme, we use a combination of 
the line search and backtracking methods, which is an efficient 
method to solve nonlinear equations. Precisely, when /j. > f^^i, 
5X may be too large. Then, 5X is reduced by a factor A, 

Xk = X^_i -H ASX, (30) 

with A chosen to minimize /(X^^i -H T^X) to maintain the itera¬ 
tion in the descending direction. 

3.2. Treatment of the singularity 

The computation of S goes through the integration of a singular 
function. Indeed, in the integrand of Eq. ( [T2| ), the function Epx) 
has a logarithmic divergence at x = 0^, see Fig[^ 

To integrate this kernel, we use the periodization method de¬ 
veloped by |Helluy et al.| ( |1998[ ), which consists of a high order 
quadrature method with a variable change. The new variable is 
a polynomial of degree k, properties of which improve on the 
order of rectangle rule. The advantage of this method is that it 
can handle singular functions. Indeed the error |£’a'(/)I on an 
interval subdivided in N subintervals goes as ~ CklN^, where 
y —» (2k - 1) if / has a logarithmic singularity at one bound 
of the integration domain (see Table [^. However the periodiza¬ 
tion method requires that the integration is done over the interval 
[0,1] and that the singularity is located at one extremum of this 
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Table 1. Benchmark of the periodization method for a logarithmic sin¬ 
gularity (from |Helluy et al.| ( |1998| l). 


f(x) 

k 

Ck 

r 

|£io(/)l 

|£2o(/)I 

ln(x) 

2 

3.22x10“ 

3.03 

0.30x10-^ 

0.37x10-5 


3 

3.34x10^ 

5.05 

0.31x10^3 

0.93x10-5 


5 

3.05x10^ 

9.09 

0.14x10-"^ 

0.42x10-^ 


10 

1.74x10''^ 

19.1 

0.48x10-5 

0.92x10-" 



Fig. 3. Map of the integrand of Jij for a radiative transition (see Eq. 
[T^. The diagonal is singular, and the quadrature in the periodization 
method produces a mesh refinement around it. The color encoding of 
the integrant is in logarithmic arbitrary units. 


interval. We thus split the integration domain into two subdo¬ 
mains, {Ty,Tv\ and [tv,t“"] , and then normalize each subdo¬ 
main to [0,1]. However, since the MPI strategy limits the num¬ 
ber of model layers handled by each process, this can lead to the 
undersampling of one of the subdomains, when Ty approaches 
one of the bounds of the total domain (i.e., the boundaries of 
the model atmosphere). We have thus tested this method under 
the same conditions as those met in our problem, without any 
resampling of the subdomains. The result is displayed in Fig. ^ 
As expected, the relative error strongly increases when one of the 
subdomains includes too few model layers (typically fewer than 
4) and becomes unacceptable for the extreme points. They will 
thus be considered ghost points in all subsequent computations. 
This does not affect the computation of the statistical equilib¬ 
rium. Because we can compute the nonlinear function F(X) with 
good accuracy, the next step to find its root is to compute the 
Jacobian matrix to solve the system through a Newton’s scheme. 

3.3. Computation of the Jacobian matrix 

The method presented here, and more generally the solving of 
nonlinear systems of equations, requires the computation of the 
Jacobian matrix. One of the advantages of the method we de¬ 
velop is that an analytical form can be derived. Here, the function 
F(X) is a set of A X P equations, N being the number of energy 
levels of the species considered, and P the number of layers in 
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Fig. 4. Benchmarking of the periodization method applied to the inte¬ 
gration with a logarithmically singular kernel {g{t) = flo + ait, K(x, t) = 
-7 - ln{\G(t) - G(jc)l) and G{x) - g(s)ds). The integration domain 
has been split into two subdomains (see text for details). 


the discretized model atmosphere. The differentiation of Fj over 
tc™ (i and k referring to energy levels, and I and m to atmosphere 
layers) leads to a set of integro-differential equations, with the 
following four cases: 


if k - i and I — m 


(31a) 


dF‘ ^ ^ / - X 


dx\ ^ dx‘ 

I ]i=i I ji=i 

ifk + i and / = m 

dF‘ dj' ! _ . 

_ ^ ~~F~r ~ ~ y\lii + -k tlcolf^kij 5 

if k - i and I + m 

dF 
dx" 


dFl ^ dj‘j 

j*i ‘ 
ifk + i and t + m 
dF' dj‘. 


(31b) 


(31c) 


(3 Id) 


Equations pia[ l and pib[ l represent local couplings. Mathe¬ 
matically, they correspond to the diagonal and block-diagonal 
parts of the Jacobian, respectively. Physically, Eq. ( [Mal i de¬ 
scribes perturbations of the statistical equilibrium of a given 
level, by modifying the population of this level, within one atmo¬ 
spheric layer. Equation pib[ ) describes the perturbations of the 
statistical equilibrium of a given level, by modifying the popula¬ 
tion of other levels, still within one atmospheric layer. Equations 
Plc| l and pid| ) correspond to off-diagonal blocks and describe 
nonlocal perturbations, i.e., the influence of the populations in 
another layer. Naturally, in these last differentiations the cou¬ 
pling is purely radiative, whereas the first and second differenti¬ 
ations include collisional processes. 

To determine the terms dx'^, one needs to differentiate 
an integral expression over the variable x™. This derivation must 


be considered a functional derivative, and be performed before 
discretization (see Appendix [B] for the demonstration). 


dJij {xk\ (tv) 
6xk (i) 


= lim 

e —»0 


Jij [Xk(Tv) + sd(Tv - s)] - Jij [Xi:(Tv)] 


(32) 


where 6(Ty - s) is the Dirac distrbution. 

In the case Ty s {I m), the derivatives for each compo¬ 
nent are 


dj-L‘-[x,,X,](Tv) 

(5x,(s) 


^ DMs)K{s, Tv) + DiI(s)Ei 


X 


.j-Out 

I ^(t)xi(t)ko(t, Tv, s)dt if k > Tv 

•Jkk 

I ({t)xi(t)Ko{t,Tv, s)dt if k < Tv 

JrV“ 

(33a) 


6jj^r[Xi,X:](Tv) 


6xjis) 


= -Dij((s)Eij 

I ^(t)xi(t)Ko(t,Tv, s)dt if k > Tv 

Jkk 

I ^(t)xi(t)Ko(t,Tv, s)dt if k < Tv 

(33b) 


X 


djfj°'nXi,Xj](Tv) 

6xj(s) 


■ as)E. 


X 


t] 

.Out 


dj;C™‘[x,,X,](Tv) 

6xj{s) 


X 


r ^ij(t)Bij(t)Lo(t, Tv, s)dt if k > Tv 
Jkk 

j ^ij(t)Bij(t)Lo(t, Tv, s)dt if k < Tv 

(33c) 

-mEij 

J ^ij(t)Bij(t)Lo(t, Tv, s)dt if k > Tv 


r 


^ij(t)Bij{t)Lo(t, Tv, s)dt if k < Tv 

(33d) 


However ^ is not differentiable at x(f) = x(tv) because of 
the singularity, therefore Eq. ( |33a| l is undetermined for s - tv- 
Eor these local terms, we thus need an approximate expression 
of the radiation field. We emphasize that we only apply this ap¬ 
proximation to the computation of the Jacobian matrix. 


3.4. Sparsity of Jacobian matrix 

In Section [XT] we argue that GMRES is a well-adapted method 
for large-scale problems if the Jacobian matrix is sparse, as it 
only requires matrix-vector products, and no matrix inversion. 
The sparsity (or density) coefficient is defined by the ratio of 
nonzero elements over the total number of elements. The Jaco¬ 
bian matrix is composed by P x E submatrix of N xN elements 
leading to NP x NP elements. The block diagonal of this matrix 
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1= l 1=2 /= 3 1-4 



1 100 200 1 100 200 


Fig. 5. Pattern of the Jacobian matrix for 4 model atmosphere layers and 
200 energy levels. 


is full because of the irreducibility of the collisional coupling 
within a given layer (no selection rules for collisions). The num¬ 
ber of nonzero elements in the block-diagonal part is f* x In 
the P(P - 1) off-diagonal submatrices, the nonzero elements are 
only due to radiative coupling between different layers, and their 
number is thus linked to the selection rules for allowed transi¬ 
tions. These rules thus lead to a sparsity of each submatrix with 
a density p^oi = transitionsThe density of the complete Ja¬ 
cobian matrix is then given by 



Fig. 7. Evolution of the Jacobian matrix sparsity according to the num¬ 
ber of atmospheric layers for different values of Pmoi- 


3.5. Jacobian block-diagonal approximation 

Assuming that tjCtv) - Trit) oc (f - ry) (first order linearization) 
when t is close to Ty, because the kernel decreases when \t - Ty| 
increases, we assume that every function at t is equal to its value 

at Ty. 

Splitting the integral over t at ry and inverting the integrals 
over f and p, [Gonzalez Garcia et al.| ( |2008| l show that the mean 
intensity can be rewritten as 

^b<^v) = , , (1 -yS^(Tv) -Mry)) (35a) 

Xj(Ty) - Xi(Ty) 

Sij(Ty) = TmaxIvPiTy), (35b) 

where y6_(Tv) and jS+(Ty) are inward and outward escape proba¬ 
bilities (defined in [Gonzalez Garcia et al.| ( |2008| l). The P(Ty) in¬ 
voke the function E 2 (x), which is not singular and differentiable 
everywhere in the domain. This enables an analytic estimate of 
the bloc diagonal terms. The pattern of the Jacobian is displayed 
inFig.|5] 


PxN^-\-P(P-l)X transitions _ 1 + Pmo\(P " 1) 


Then, as PmoiiN) decreases rapidly with N (e.g., PhjO ~ 10 



Fig. 6. Example of the evolution of the off-diagonal, block-matrix spar¬ 
sity (Pnioi) according to the number of considered levels. This example 
is computed for the 411 levels of the ortho-H20 molecule. 


for A = 411, see Figs.|^and|^ and p is inversely proportional to 
P, the larger the system, the sparser the Jacobian matrix. The 
GMRES scheme is thus well adapted and asymptotically ap¬ 
proaches the NP complexity. 


3.6. The numerical code MOrad 

The above strategy has been implemented in a numerical code 
named MOrad written in Fortran 90. It takes as inputs: 

1. A model atmosphere (temperature, atomic, molecular, and 
electronic densities etc. as a function of geometrical thick¬ 
ness) 

2. Spectroscopic data (level energies, statistical weights, radia¬ 
tive rates) 

3. Collisional rates 

Then it computes the continuous opacities and uses an ini¬ 
tial guess to calculate the function ( |22) and its Jacobian ma¬ 
trix. Choosing the initial guess Xq is an important operation. In 
MOrad we tested the use of an LTE solution (Boltzmann distri¬ 
bution), and a 0 K solution (the molecules in the ground level). 
We find that the convergence rate is better with the "0 K" initial 
guess. 

The singular integrations needed to compute S are per¬ 
formed with the periodization method. The Green’s function 
needed to compute S is stored at each point. The required space 
scales as the number of layers times the number of frequency 
points (in a parallel code; otherwise it would scale with the 
square of the number of layers). The evaluation of the function 
and the Jacobian are parallelized with MPI. At first glance, two 
parallelization strategies may be considered: either a paralleliza¬ 
tion over frequencies or over atmosphere layers. However, [van[ 
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|Noort et ar] ( |2002[ ) show that the most efficient parallelization to 
reduce communications is spatial parallelization because of the 
strong local coupling between energy levels. The populations of 
all energy levels within each atmosphere layer should thus be 
computed by the same processor, with several adjacent layers 
possibly treated by the same processor. In MOrad we adopt a 
subdomain decomposition where each processor computes the 
physical values needed in one given layer and solves the local 
statistical equilibrium. Integrations are performed with a scal¬ 
able algorithm, which minimizes global communications to pre¬ 
serve a good scalability. 

The code has been interfaced with a nonlinear solver of 
the parallel and scalable library PETSc (see, e.g., [Balay et al.| 


|2013b|a[|1997| l. This solver takes as input the distributed func¬ 
tion vector and the Jacobian matrix and returns the root of the 
function after preconditioning. 

The main advantages of this library are its performance and 
its scalability. Moreover, PETSc provides access to many dif¬ 
ferent nonlinear methods and preconditioning, which allows us 
a versatile choice of the best method to be adapted for a given 
physical problem. Other efficient libraries, such as the parallel 
lO library HDE5, are interfaced with MOrad to preserve the per¬ 
formances for large-scale problems and the portability of output 
data. Euture planned developments include an MPI/OpenMP hy¬ 
bridization (openMP parallelization over frequencies). 


4. Discussion 

4.1. Example of application 

To illustrate the applicability of MOrad, we compute the nonLTE 
departure coefficients of water molecules in a MARCS ( [Gustafs-] 
|son et aL]|2008| l red supergiant model atmosphere. The model 
considered here has an effective temperature of 3500 K, a logg = 
0, a solar metallicity, and a microturbulence of 2 km s ^ These 
parameters correspond to typical red supergiant stars. We con¬ 
sider more than 800 rovibrational levels, leading to more than 
330 000 transitions and 15 000 lines (see the Grotrian diagram 
of ortho H 2 O in Eig.[^. The energy levels and the radiative co¬ 
efficients are taken from [Barber et al.| ( |2006|l. The collisional 
rates for H 2 O-H 2 and H20-e“ are taken from jEaure & Josselinj 

( I 2 OO 8 I 1 . 

5,000 


4,000 

E 

— 3,000 

2,000 

> 

J 

1,000 


0 



Fig. 8. Grotrian diagram of the ortho-H20 molecule showing the ra¬ 
diative transitions taken into account. The parameter J„ is the rotational 
angular momentum of H 2 O. Colors indicate different vibrational states. 

The code MOrad was launched on 47 processors using a 
GMRES Newton-Krylov method and a line-search global con¬ 
vergence strategy. We conducted preconditioning with the AMS 
method (Block Additive Schwarz method), where each subblock 
is preconditioned with an approximate iterative LU factorization 


method. We obtain the root X* of the function (||F(X*)||oo < 
10 '^) after three nonlinear iterations and a total time of code 
execution of less than ~lh on the Alarik LUNARC system (see 
|Lunarc|2013| l. The departure coefficients are presented in Eig.|^ 
The convergence rate is Q-quadratic, though a quadratic conver¬ 
gence rate would be expected. This may be because of an over 
reconditioning or a problem that is too simple(close to linear), 
but we obtain the same Q-quadratic convergence for a purely ra¬ 
diative problem. 



10 > 10 ’ 10 ‘ 10“ 

Optical Depth Ty [vy = 500 nm] 


Fig. 9. LTE departure coefficients of ortho-H20 in a red supergiant 
model atmosphere. 

A detailed analysis of the results is out of the scope of this 
paper and is devoted to a forthcoming paper; a preliminary anal¬ 
ysis was done in [Lambert et al.j ( [2013[ ). In particular a detailed 
discussion of the use of the super-level approximation, which 
seems appropriate for the vibrationally excited states, will be 
presented. In short, the non-LTE calculations lead to stronger 
lines, especially rotational lines in the fundamental vibrational 
state around 2 pm, compared to LTE computation. 

4.2. Performance and complementarity of the method 

The method presented here is naturally accurate as it is equiva¬ 
lent to finding the root of a function and thus avoids any false 
convergence (assuming the solution is unique^ Eurthermore, 
we emphasize that because the integration is analytic, the method 
leads naturally to null residue. Concerning the performance of 
the method, the performances of the method are difficult to as¬ 
sess. The convergence rate depends upon the choice of both the 
solver and the preconditioning. These choices are not unique 
and may vary according to the considered problem (species and 
model atmosphere). The memory demands scales with the di¬ 
mensionality of the considered problem. 

The computation of the function and the Jacobian matrix 
leads to a scheme with an asymptotic complexity in OiNb^sq x 
{N X P)^). The complexity of the resolution scheme is more 
difficult to evaluate because the preconditioned Newton-Krylov 
method has a cost in Afreq x x (Nfreq being the number 
of frequency points sampling the line profile, N the number of 
energy levels, and P the number of layers in the discretized at¬ 
mosphere), but invokes another term that evolves as the square 

^ The solution of the Eddington problem (the two-level atom in the 
isothermal optically thick medium) is recovered at the 1st iteration with 
preconditioning, and after five iterations without preconditioning (rel¬ 
ative error < 0.005 for e = 10“^ and = 1000). This is expected 
because the inversion of the matrix by the Krylov subspace method is 
well known to work in this kind of situation. 
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of the iteration index. This term is partially eliminated by restart¬ 
ing GMRES and is not really important in practice because the 
number of iterations is small after appropriate preconditioning. 

Furthermore, the performance is better than (9(n^) and prob¬ 
ably close to 0(n). Indeed, for a sparse matrix, GMRES evolves 
in 0{n). Because the only nonzero elements of the off-diagonal 
blocks of the Jacobian matrix are due to radiative coupling be¬ 
tween layers, the sparsity of J evolves as the ratio of the num¬ 
ber of radiative transitions over N^. Taking a complexity in 
Nbfreq X X is then an superior limit, which underevaluates 
the method scalability in real cases. 

As a comparison, the computation of the equivalent function 
of F(X) in MULTI has a complexity evolving in Ax/*because of 
the Scharmer operator, which performs no numerical integration 
to evaluate diffusion integrals but only one quadrature point. The 
system is solved with a Newton-Raphson scheme, which has a 
cost in X when using the full matrix or x P when us¬ 
ing the local operator (block-diagonal matrix; M. Carlsson, priv. 
comm.). The theoretical asymptotic speed-up for MOrad com¬ 
pared to MULTI are summarized in Table |4.2| 


fMULTl/fMOrad 

Full Operator OAB (diagonal) 

p = 1 N N/P 

p = 0 PN^ 


Table 2. Theoretical asymptotic speed-up for two extreme cases in 
terms of p = sparsity(J)(0 < p < 1) for MOrad, compared to MULTI 
with either full or local operator (OAB stands for Olson, Auer, & Buch- 
ler Operator). 


From Table 4.2 one can see that our method is complemen¬ 
tary to existing methods, such as direct linearization. Indeed, in 
some cases, with a large number of levels and a small number of 
layers, the expected time of computation is shorter. 


5. Conclusion 

We have developed a new method to solve the RTF in non-LTE 
conditions. Rather than solving the RTF itself, we solve the sta¬ 
tistical equilibrium equation system, which includes an analyti¬ 
cal formal solution of the RTE. One of the main advantages of 
this choice is that the coefficients of the associated Jacobian ma¬ 
trix can be computed exactly. With this approach, a nonlinear 
solver can be used. In other words, this method avoids any lin¬ 
earization or stationary iterative schemes, which may converge 
very slowly or even falsely in some extreme cases (a very low 
non-LTE parameter e and/or highly discretized media). We chose 
a method based on the GMRES nonlinear method. It has been 
implemented in a code that uses the PETSc library. We pay spe¬ 
cial attention to the proper treatment of the singularity that arises 
in the integration of the mean radiation field. 

Furthermore, this method is parallelized using MPI, which 
makes it applicable to large-scale systems such as molecules and 
atoms in stellar atmospheres or interstellar medium, as it has an 
excellent scalability with the number of energy levels. Indeed, 
our code has been successfully applied to the modeling of water 
in the atmospheres of red supergiants. The results will be pre¬ 
sented in a forthcoming paper. Further developments of the code 
include hybrid parallelization, with openMP parallelization over 
frequencies. We then plan to make the code public. 
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Appendix A: Formulation 

One starts with the classical formulation of the radiative transfer equation in a plane parallel geometry, 
dly ( 


Cont , -^Line 
Xv "•"/V 


)/v 


+ 


(A.l) 


The monochromatic extinction and emissivity coefficients are expressed with an explicit dependence on levels populations as fol¬ 
lows: 


.Line _ ... n ^ \ 

Xv ~ Xul ~ Yv A 
47r 

„Line _ iid^ _ jul . 

~ Yy ^ul ~ Yy a ^ul^U. 

4;r 


Tv'“ = 

4;r 

„Line _ a r 


(A.2a) 

(A.2b) 


By inserting relative populations / = tiiln^ .where is the total density of the species labeled by k, one gets 

(A.3a) 
(A.3b) 

One brings up the new variable x, = filgi, and using the relation Buigu - Biugi, where g, is the statistical weight of the level i, yields 


Line k /n rt k n / \ 

Ay = -T—n guiBiu - Bui—) = 0y -;—n guBuiixi - x„) 

47r gu gu 47r 


Line lul k a fu 

dv guAui — 

47r gu 


K'^^d^guAulXu. 

4;r 


Including the formulation of;^fJ;“® and rj^ in the radiative transfer equation, one obtains 


B 


dly 

ds 


Ay°"‘ + K'’^n^guBui(xi - x„) ly + + cPl^'^n'^guAuiXu. 


An 


An 


(A.4a) 

(A.4b) 

(A.5) 


In order to create a transition independent scale, the previous equation is divided by an arbitrarily chosen continuum extinction 
coefficient ;^fy, here;^fy = Xv= 500 nm- Moreover we express it per hydrogen atom to bring out the ratio n^/n^ for numerical reasons. 
This normalization is symbolized by a tilde (~) for other variables. 


dly 

1 - 

Xvds 


hvui n'^ 1 


A H ~ ^ulSu i.^1 -^m) I , jj _ ^ulSu^t 


Xv 4;r Xv 


hVui 1 


y-Cont „Hr;>Cont 
juI , tv Xv 


An Xv 


n^^Xv Xv 


.Cont 


(A.6) 


With ^y = 1^, Eij = i^Bijgi, ( = ^ Oij = i^Aijgi, Sy = and dry = -Xvds. Assuming that the continuum is formed ir 
LTE = By), one obtains 

dL 


= [^V + Eul^4>y(Xl - Xu)] ly - DulXu((py - ^yBy. 


(A.7) 


Appendix B: Jacobian matrix coefficients 

We start with the analytic formulation of the mean intensity field for the line component (similar demonstration for contunium), i.e.. 


Jij[Xi,Xj]{Ty) = D, 


•f 


^(t)Xi(t)K[Xi, Xj](t, Ty)dt. 


Because the functional derivative of a functional product is 

(5(E[y]G[y]) /E[y] /G[y] 

-j-TT - = G[y]—+ E[y]——, 

6y{x) 6y(x) 6y(x) 

we functionally derive the product Xi{t)K[xi, Xj\ over .r,(Tv = s), which yields, 

djij[Xi,Xj](Ty) 


(B.l) 


(B.2) 


6xi(s) 


= A 


r 


6xi(t) 

Xj\(t, Ty)dt + Dij 


r 


^(OXi(t) -——- dt. 

OXi(s) 


(B.3) 
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Equation ( |B.3| l can be rewritten as a sum of two integrals over f, Ii and /a. 
The first term /i simplifies 


h = Di 


'f 

*JT,, 


6Xi(t) 

^(0 ^ ^^ K[xi,xj](t,Tv)dt - Dij 


6xi(s) 


f 


^{f)6{t - s)K{xi, Xj\{t,T\/)dt 


= Dij((s)K(s,Ty) 


(B.4) 


This term is an indefinite form if i = xy, because the kernel K{x,y) is singular atx—y. 

We consider here the case s xy (the case i = ry is detailed in Section 3.4). 

The second term I 2 requires the computation of the functional derivative of the kernel, which is 

K[Xi, Xj](t, Tv) = - J tpy(Xy)(py{t)Ei Xj]{t, Ty)) dv. 

With the chain rule theorem of functional derivatives, we get 

mm) ^ df{F\y])6F\y] 

6y(x) dF\y] 6y(x)' 

We then derive the kernel K, finding 


(B.5) 


(B.6) 


6K[Xi,Xj]{t,Xy) 1 


6xi(s) 


----- r 

2 Jo 


4>y(Xy)^y(t) X 


•*)](?, Ty)) 


SXiis) 


dv 


1 

^~2 J X Fo[Ax1°'-[Xi,Xj\{t,Xy)) 


6Axl°'^[Xi,Xj\{t,Xy) 

6xi{s) 


dv. 


(B.7) 


To simplify the calculus, we rewrite as 

^Tv 


AT^°‘[x,-,xd(f,Tv) = 


-£ 


^y(s') + £';/(s')0v(s') ds' 

@(t, Ty, /) [^v(s') + Fij^(s')(py(s') [.r/s') - X,■(/)]] ds'. 


(B.8) 


where the distribution 0(f, Ty, i) = -2(6(t - xy) - O.5)(0(s -t)- 9{s - Ty)), and 9{x) is the Heaviside function. This distribution can 
be rewritten. 


0(f,Ty,s) = ( . 

if S > Ty 

We functionally derive Eq. ( |B.8| ), yielding 
5ATj°‘[Xi,X;](f,Tv) 


(B.9) 


- I &{t,xy,s')Fij^{s')(l>y(s')^p-^ds' = - f &{t,xy,s')Fim)(py{s')6{s - s')ds' 
6xi(s) Jyln 6xi(s) Jyin 

= -Fij^(s)&(t,Xy, s)4>y(s). 


(B.IO) 


We insert this expression in Eq. (B.7 1 and we define the modified kernel Ko(t, xy, i), 

6K[Xi, Xj](t, Xy) 


6xi(s) 


1 

Fij((s)@{t,Xy, s)- J c)>y(xy)(py(t)M'^)Eo(^T'^°Jxi,Xj](t,Xy)jdv 


F,j{(s)0(t,Xy, s)Ko[Xi,Xj](t,Xy, s). 


(B.ll) 


The integral I 2 from equation ( |B.3| l can be rewritten as 

h = D,j 


'f 


6K[Xi, X;](f, Ty) 
^(t)Xi(t) - ; ! dt 


DijFms) 


r 


6xi(s) 

&(t, Xy, s)^(t)Xi(t)Eo[Xi, Xj](t, Xy, s)dt, 

which, according to the 0 properties, becomes 

^(t)Xi(t)Ko{t,Xy, s)dt if i > Ty 
^{f)Xi(t)Ko{t, Xy, s)dt if S < Ty 
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h = DijEij({s) X < 


r 

r 

kJt,, 
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Finally, we get the derivative of the mean radiation field, 


6Jij[Xi,Xj](TY) 

6xi(s) 




K(s, Tv) + Ejj X < 


r 

X 

kJt., 


((t)xi(t)KQ(t,Tv, s)dt 
^(t)Xi(t)Ko(t,Ty/, s)dt 


\ 

if S > Tv 
if i < Tv 


Appendix C: Definition of the variabies 


(B.14) 


Table C.l. Physical quantities 


Symbol 

Meaning 

units (cgs) 


Direction cosine ip. - cos(0)) 



T 

Temperature 

K 



Total density of species k 

cm“^ 


Hi 

Density in an energy level labeled by i 

cm“^ 


f 

Relative density {f, - riijn'^) 



Xj 

Weighted relative density {xi - filg,) 



8i 

Statistical weight of level i 



V 

Frequency 

Hz 


ly 

Specific intensity 

ergs- 

' cm ^ sr 'Hz ' 

Jy 

Mean intensity 

ergs- 

' cm-^ Hz-' 

Ty 

Optical depth 



Tv 

Optical depth at 500 nm 



Sy 

Source function 

ergs- 

' cm-^ sr-' Hz-' 

By 

Planck function 

ergs- 

' cm-^ sr-' Hz-' 

Xy 

Monochromatic extinction 

cm-' 


Tjy 

Monochromatic emissivity 

ergs- 

' cm-^ sr-' Hz-' 

Py 

Escape probability 



Aut 

Einstein coefficient for spontaneous de-excitation 



Bu 

Einstein coefficient for stimulated (de)excitation 

erg-' 

cm^ Hz 

Cij 

Rate for collisional (de)excitation 

cm^ s 

-1 

4>y 

Normalized line profile of the transition i —» j 

Hz ' 



Table C.2. Functions, operators and distributions 


Eunctions 

Name 

Meaning 

E„ix) 

n''^ Exponential integral 

/ e 

En - 1 - dt 

1, f" 

By(T) 

Planck function 

2hv^ 1 

By(T) - ^2 ^l,y/k^T _ 1 

6(x) 

Dirac distribution 

« = \o 



XX, 6{x)dx - 1 

n(x) 

Gate distribution 

n(;r) = H(x + 1/2) - H(x - 111) 


Generalized gate distribution (boxcar) 

'^a,bix) - n((x -ib + a)l2ib - a))) 

e{x) 

Heaviside distribution 

|l if ;c > 0. 


Green’s function 

-C^ix, s) - Six - s) 

1 

Integration of / weighted by a line profile 

7 = r fy<f>ydv 

\J —oo 
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